Pairing in population imbalanced Fermion systems 
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We use Quantum Monte Carlo (QMC) simulations to study the pairing mechanism in a one- 
dimensional fermionic system governed by the Hubbard model with attractive contact interaction 
and with imbalance between the two spin populations. This is done for the uniform system and 
also for the system confined in a harmonic trap to compare with experiments on confined ultra- 
cold atoms. In the uniform case we determine the phase diagram in the polarization-temperature 
plane and find that the "Fulde-Ferrell-Larkin-Ovchinnikov" (FFLO) phase is robust and persists 
to higher temperature for higher polarization. In the confined case, we also find that the FFLO 
phase is stabilized by higher polarization and that it is within the range of detection of experiments 
currently underway. 



The best understood mechanism for pair formation of 
fermions is the BCS mechanism 1 - where two fermions with 
opposite spin and equal but opposite momenta form a 
pair with zero center-of-mass momentum. Shortly after 
the development of the BCS theory of superconductivity, 
the question of pair formation in polarized superconduct- 
ing systems, i.e. where the populations of the two spin 
states are imbalanced, was addressed independently by 
Fulde and Ferrel 2 - (FF), Larkin and Ovchinnhkm* 3 - (LO) 
and Sarma 4 -. Initially, the question was motivated by 
interest in the nature of superconductivity in the pres- 
ence of a magnetic field but since then other instances 
where such a mechanism intervenes have become of in- 
terest. For example, in the interior of supermassive 
neutron stars, quarks of various colors may form pairs 
which are not colorless thus leading to what is known as 
"color superluid"— . Another situation of major current 
experimental interest is in systems of confined ultra-cold 
fermionic atoms such as 6 Li or 40 K. Such experiments 
have now reported the presence of pairing in the case of 
unequal populations^ in three-dimensional cigar shaped 
traps and in one dimensional traps^. However, the pre- 
cise nature of the pairing has not yet been elucidated 
experimentally. 

On the theoretical side, many methods have been used 
ranging from mean field^ to effective Lagrangian 1 ^ to 
Bethe ansataii. The two competing mechanisms for 
pair formation in the population imbalanced case are 
the FFLO and the Sarma mechanisms. In the former, 
the bosonic pairs form with non-zero center of mass mo- 
mentum equal to the difference in the Fermi momenta, 
\kpi — kp-2\, of the two populations. This leads to the for- 
mation of a standing wave in the order parameter whose 
wave vector is \kpi — &F2|- With the Sarma mechanism, 
majority fermions whose momenta equal the Fermi mo- 
mentum of the minority, are promoted to higher momen- 
tum levels thus forming a breach in the Fermi distribution 
of the majority population. This breach allows majority 
fermions with momentum equal to the Fermi momentum 
of the minority to pair up with minority fermions form- 
ing pairs with zero center of mass momentum. Exten- 
sive numerical work on the one dimensional system us- 



ing Quantum Monte Carlo (QMC)^^, and the Density 
Matrix Renormalization Group (DMRG) 1 ^) has demon- 
strated that, in the ground state, population imbalance 
leads to a robust FFLO phase over a very wide range of 
polarization and interaction strengths. The Sarma phase 
was not detected in these numerical works. 

The stability of the FFLO phase at finite temperatures 
has been adressed with mean field calculations 1 ^ which 
can be unreliable in low dimension where quantum fluc- 
tuations are large. Recently we addressed this question 1 ^ 
using exact quantum Monte Carlo (QMC) simulation 
which is the focus of this presentation. 

In order to study the pairing mechanism of fermions 
in an optical lattice, we consider the one-dimensional 
fermionic Hubbard Hamiltonian, 
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where c\ a and c ia are fermion creation and annihilation 
operators on lattice site i satisfying the usual anticom- 



mutation relation, {c 



The fermionic 



species arc labeled by a = 1,2 and = c\ a c ia is 
the corresponding number operator. The energy scale 
is set by taking the hopping parameter 4=1. The con- 
tact interaction strength U is negative since we are in- 
terested in pair formation in the attractive model. The 
last term describes the confining harmonic trap which is 
centered at the midpoint, L/2, of the i-sitc lattice. We 
take periodic boundary conditions. Our QMC results 
were obtained using the Determinant QMC algorithmic 
(DQMC) and the Stochastic Green Function (SGF) 
technique 1 ^. Specifically, we used the DQMC algorithm 
to determine the phase diagram in the polarization- 
temperature (P, T) plane because of its good convergence 
properties for large systems. This algorithm functions 
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in the grand canonical ensemble where the chemical po- 
tential for each species is tuned to obtain the desired 
polarization. When the populations are imbalanced, this 
algorithm suffers from the sign problem even with the at- 
tractive interaction, as is the case here. However, in our 
simulations, the average sign never went below around 
0.3 at the highest polarizations and lowest temperatures. 
This allowed us to study the system under rather ex- 
treme conditions. The SGF algorithm was used for all 
other simulations including all simulations of the confined 
system because this algorithm functions in the canonical 
ensemble, where the populations are fixed, which corre- 
sponds to the experimental situation. In addition, we 
have verified^ that for the system sizes and fillings we 
simulated, the grand canonical (DQMC) and canonical 
(SGF) ensembles give the same results. Typical simula- 
tions, with DQMC or SGF, took from five to seven days 
each on a 3GHz processor. The simulations were per- 
formed locally on our cluster with 102 cores. 




FIG. 1: Top: The momentum distributions for the majority, 
minority and pair populations in a system of L = 32 sites, 
majority population Ni = 9, minority N2 = 7, /3 = 64. The 
pair momentum distribution exhibits a distinct maximum at 
nonzero momentum. Bottom: The pair momentum distribu- 
tion for several polarizations P. The inset shows that the 
position of the FFLO peak is equal to \kpi — &F2|. 

Using these algorithms, we calculate the real space 
Green functions of each species, G a , and the pair green 
function, G n 



r pair i 



G 



pair (0 

A, 



^3+1 <r C j a' 
Cj2Cjl, 



(2) 
(3) 
(4) 



where Aj destroys a pair on site j. The Fourier transform 
of G a (l) yields the momentum distributions n a (k) and 
the transform of G pa ; r (7) leads to the pair momentum 
distribution, n pa i r (fc), a central quantity in this work. 
In the non-interacting limit, the Fermi momentum of a 
population is given by kp a = Na 2 ~ 1 T" ' wncrc L ls the 
number of sites and N a the number of particles. The 
Fermi energy is given in the uniform case by ep = tkp = 
Tp where Tp is the Fermi temperature. The polarization 
is defined as P = {N\ - N 2 )/N where N = Nj_ + 7V 2 is 
the total number of particles. 



In Fig[T] we show results for the uniform systei 
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Vt = 0. The top panels show the momentum distri- 
butions for the majority and minority populations and 
for the pairs at U = —it, — lOi. We see that the Fermi 
"surfaces" for the two populations are sharply defined 
and that the peak in the pair momentum distributions is 
at nonzero momentum, k pca ^: This is a clear signal for 
FFLO pairing. It is interesting to note the deformation 
of the momentum distribution of the majority popula- 
tion: A bump develops in n\(k) for k > kp2- As the 
Fermi distributions of the minority and majority start 
to match up more closely as U becomes more attrac- 
tive, the excess unpaired majority particles populate the 
states with momenta larger than the minority Fermi mo- 
mentum. This give the bump in n\(k > Kp2). The 
lower panel of FigQ] shows n pa i T (k) for several polariza- 
tions. It is clearly seen that as P increases so does A; pea k, 
the momentum at which n pa i r (fc) peaks. In the inset we 
display fc pca k versus |fcpi — &F2I which shows that in fact 
&pcak = I&F1 — kp2\ as predicted by the FFLO scenario. 
The nonzero value of the pair center of mass momen- 
tum means that its Fourier transform G pa i r (7) (Eq. [3])) 
oscillates^ with wavelength A = 27r/|fc pea k|. This means 
that in the FFLO phase the system is not homogeneous, 
it consists of pair-rich regions and regions depleted in 
pairs but rich in the excess unpaired majority popula- 
tion. 

Such calculations were done for several values of U 
using DMRG— and QMCi2 and lead to the conclusion 
that, in the ground state, the polarized system is always 
in the FFLO phase and the Sarma mechanism does not 
intervene. The question then arises as to the effect of 
finite temperature: Will FFLO survive at T 7^ and 
how robust is this phase? 

To this end, we studied, using QMC simulations, what 
happens to the FFLO peak as the temperature is in- 
creased. The inset in Fig. [5] shows this peak for the case 
Ni = 13 and 7V 2 = 7 (P = 0.3). We sec that as T in- 
creases, the FFLO peak gets lower and eventually disap- 
pears at p = 1 /T « 5 which signals the destruction of the 
FFLO phase by thermal fluctuations. But how is FFLO 
destroyed? The two possibilities are (a) the pairs are bro- 
ken or (b) the pairs are still formed but thermal fluctua- 
tions make the system homogeneous resulting in a peak 
at k = for n pea k(fc)- To destroy the pairs, the thermal 
energy should be of the order of the pair binding energy, 
i.e. T ~ \U\. However, we see from Fig. [2] that FFLO 



O.N 



0.6 



0.4 



0.2 




Q-0|3=32 
h-hP=16 

^P=10 
<J-< P=5 

0-0 P=3 



I 


i 1 i 






i 



0.1 



0.2 



0.3 



0.4 



0.5 













O-O P=0.03 






as P=0.05 






0-0 P=0.08 


- 




A-A P=0.1 1 


- 




<-< P=0.15 






V-VP=0.18 






L> — t> V=\J.eLtL 






-1— <r P=0.26 






x-x P=0.3 






*-» P=0.35 






•-• P=0.39 






■ « P=0.44 




b 


♦-♦ P=0.5 




P=0.56 



FIG. 2: The phase diagram at quarter filling, N = L/2, with 
U = — 3.5t and L — 50. The Fermi temperature at this filling 
is Tp = 0.5t. We see that the FFLO phase is robust and 
persists over a wide range of P and T. Increasing P stabilizes 
FFLO up to higher T. PPP stands for Polarized Paired Phase, 
see text for a discussion. The inset shows the behavior of the 
FFLO peak as T is increased for the case iVi = 13, N2 = 7, 
i.e. P = 0.3. We see that the FFLO peak goes down as T 
increases and vanishes at f3 — 1/T ~ 5. 



is always destroyed at T <C \U\. We conclude, there- 
fore, that thermal fluctuations destroy FFLO by making 
the system homogeneous not by destroying the pairs; we 
denote this phase by Polarized Paired Phase (PPP). In 
this way, we determine the phase diagram of the system, 
shown in FigE] at quarter filling N = Nx+ N 2 = L/2. 
The figure shows clearly that the FFLO phase is very ro- 
bust extending over a wide range of P and T . Increasing 
P stabilizes FFLO up to higher T while at low P even a 
small increase in T destroys it. Consequently, for small 
P one needs to simulate the system at exceedingly low 
T to see FFLO. This increases the simulation time and 
sets a limit on the lowest practical P. The dotted line in 
Fig. [2] schematizes the phase boundary at very low P. At 
the other extreme, P — > 1, the system is primarily made 
of one population with very few minority particles. This 
makes the FFLO signal very difficult to see. The open 
symbols in Fig. [2] denote the highest P we were able to 
examine, the system up to those values is still FFLO. 

Continuing earlier work in higher dimension^, the 
Rice group^ recently reported on experiments in arrays of 
one dimensional tubes of confined Fermionic atoms ( 6 Li) 
with imbalanccd populations. Along the tube, the atoms 
were confined with a trap frequency oj z = 2tt x 200Hz; in 
the central tube, the total number of atoms at zero po- 
larization was approximately 250 and the temperature 
was estimated at T/Tp w 0.1 where the Fermi tem- 
perature Tp is obtained from the Fermi energy ep = 
(1/2 + Huj)N/2. The pair binding energy, e = h 2 /mam 
(where cl\d is the effective one-dimensional scattering 
length), was estimated to be e/ep w 5.3. 



FIG. 3: The pair momentum, 



r (fc), versus k for several 



polarizations of the confined system, Vt = 0.0007t, L = 120, 
U = —Wt and /3 = 64. The majority population is Ni = 39 
and the minority, N2, is varied to tune P. Inset: Profile of 
the populations difference for the smallest four polarizations. 
See the text for a discussion of the oscillations. 



We now present QMC results for the fermionic Hub- 
bard model, Eq.([l|), in the presence of the confining trap. 
To make contact with the experiment we introduce the 
trapping potential in Eq.([T|) Vt = 0.0007£ which corre- 
sponds to huj z = 2\JWt- The total number of particles in 
our simulations for balanced populations is 78, to be com- 
pared with 250 in the experiment. We performed our sim- 
ulations in the temperature range 0.008 < T/Tp < 0.25 
which includes the temperature at which the experiments 
were performed, T/Tp = 0.1. In addition, to place our 
system in the same coupling parameter regime as the ex- 
periments, we present our results for a coupling strength 
of U = — Wt. U is the "pair binding energy" and the 
value we have chosen gives \U\/ep = 4.8, close to the 
experimental value. 

In FigJ3]we show, as we did in the uniform case FigJlJ 



the pair momentum distribution, 



r (fc) for several P 



values. We see that fc pea k 7^ and that its value in- 
creases with P as predicted by the FFLO picture. These 
results were obtained at /3 = 64 and essentially represent 
the ground state behavior of the system. Therefore, the 
presence of the trap does not change the nature of the 
phase in the ground state, it remains FFLO when the 
system is polarized and is robust. The inset in Fig. [3] 
shows the difference between the majority and minority 
density profiles for the four lowest P values we examined. 
We see that the difference is oscillatory; the wavelength 
of the oscillations is in fact given by A = 27r/|fcpi — fc F2 | |. 
For example, we see that for the smallest P we show, one 
wavelength fits in the system and for the P = 0.11 case 
four wavelengths fit. This is a nice visual confirmation 
that the FFLO state is not homogeneous. 

The effect of finite temperature is examined, as before, 
by studying the behavior of the FFLO peak as j3 = 1 /T 
is decreased (increased). We show in Fig0]two P cases. 
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FIG. 4: Density profiles for P = 0.05 at T/T F = 0.01 (a) 
and T/Tf = 0.1 (b). The pair momentum distribution for 
three T values (c) shows that the FFLO peaks disappears for 
T/Tf — 0.018 which is much lower than what is feasible ex- 
perimentally. The high polarization case, P = 0.56 is shown 
in (d) and (e) for T/T F = 0.02 and 0.11. The pair momen- 
tum distribution (f) shows that for P — 0.56 the FFLO peak 
survives even for T/Tf > 0.1, the experimental value. 



On the left, (a) and (b) show how the density profiles, 
i.e. the local density in the trap, for P = 0.05 change as 
the system is heated. We see that the profiles get more 
rounded as T increases. The pair momentum distribu- 
tion (c) shows that the FFLO peak disappears by the 
time T = 0.018Tp. This temperature is very low and 
is not accessible experimentally. However, the high po- 
larization case, P = 0.56, shown in (d) and (e) behaves 



differently. Here too, the profiles get rounded as T in- 
creases, but we see in (f) that the FFLO peak survives 
for T > O.lTp, the experimental value. This means that 
with presently attainable experimental temperatures, the 
FFLO phase may be observed. This increased stability 
of FFLO with increased P is consistent with the phase 
diagram we found for the uniform case, Figj2j 

In this paper we have examined the pairing mecha- 
nism in fermionic systems with imbalanced populations 
both in the absence and presence of a confining trap 
which breaks translational invariance. We showed that 
the dominant pairing mechanism in the ground state 
and also at finite temperatures is FFLO where the pairs 
form with a nonzero center of mass momentum. This 
is revealed clearly by a peak at nonzero momentum, 
^pcak = I^Fi — in the pair momentum distribution. 
The behavior of this peak is studied as a function of the 
temperature and also the polarization, P. We showed 
that increasing P stabilizes FFLO up to higher temper- 
atures and, in the confined case relevant to experiments 
on ultra-cold fermionic atoms, places this phase within 
reach. 

Acknowledgments This work was supported by: an 
ARO Award W911NF0710576 with funds from the 
DARPA OLE Program; by the CNRS-UC Davis 
EPOCAL joint research grant; by NSF grant OISE- 
0952300; by the France-Singapore Merlion program 
(PHC Egide, SpinCold 2.02.07 and FermiCold 2.01.09) 
and the CNRS PICS 4159 (France). Centre for Quantum 
Technologies is a Research Centre of Excellence funded 
by the Ministry of Education and the National Research 
Foundation of Singapore. 



1 J. Bardeen, L.N. Cooper and J.R. Schrieffer, 
Phys. Rev. 108, 1175 (1957). 

2 P. Fulde and A. Ferrell, Phys. Rev. 135, A550 (1964). 

3 A. Larkin and Y.N. Ovchinnikov, Zh. Eksp. Teor. Fiz. 47, 
1136 (1964) [Sov. Phys. JETP 20, 762 (1965)]. 

4 G. Sarma, Phys. Chem. Solids 24, 1029 (1963); S. Takada 
and T. Izuyama, Prog. Theor. Phys. 41, 635 (1969). 

5 V.L. Ginzburg and D.A. Kirzhnits, Soviet Phys. JETP 20, 
1346 (1965); Nature 220, 148 (1968); R. Casalbuoni and 
G. Nardulli, Rev. Mod. Phys 76, 263 (2004). 

6 M.W. Zwierlein it et al., Science 311, 492 (2006); 
M.W. Zwierlein et al, Nature 422, 54 (2006); Y. Shin et al, 
Phys. Rev. Lett. 97, 030401 (2006); Y. Shin, C.H. Schunck 
et al, Nature 451, 689 (2008). 

7 G.B. Partridge et al, Science 311, 503 (2006); G.B. Par- 
tridge et al., Phys. Rev. Lett. 9 7, 190407 (2006). 

8 Y. Liao et al. larXiv:0912. 0092K 4 [physics, atom-ph]. 

9 P. Castorina et al, Phys. Rev. A72, 025601 (2005); 
D.E. Sheehy and L. Radzihovsky, Phys. Rev. Lett. 96, 
060401 (2006); J. Kinnunen, L.M. Jensen and P. Torma, 
Phys. Rev. Lett. 96, 110403 (2006); K. Machida, 
T. Mizushima and M. Ichioka, Phys. Rev. Lett. 97, 120407 



(2006); M.M. Parish et al, Nature Phys. 3, 124 (2007); 
H. Hu, X.-J. Liu and P.T. Drummond, Phys. Rev. Lett. 98, 
070403 (2007); T. Koponen et al, New J. Phys. 8, 179 
(2006); X.-J. Liu, H. Hu, P.D. Drummond, Phys. Rev. 
A76, 043605 (2007). 

D.T. Son and M.A. Stephanov, Phys. Rev. A74, 013614 
(2006). 

1 G. Orso, Phys. Rev. Lett. 98, 070402 (2007). 

2 G.G. Batrouni, M.H. Huntley, V.G. Rousseau and R.T. 
Scalettar, Phys. Rev. Lett. 100, 116405 (2008). 

3 M. Casula, D.M. Ceperley, and E.J. Mueller, Phys. Rev. 
A78, 033607 (2008). 

4 A. Feiguin and F. Heidrich-Meisner, Phys. Rev. B76, 
220508(R) (2007); A. Liischer, R.M. Noack, and 
A.M. Lauchli, Phys. Rev. A78, 013637 (2008); M. Rizzi 
et al, Phys. Rev. B77, 245105 (2008); M. Tezuka and 
M. Ueda, Phys. Rev. Lett. 100, 110403 (2008). 

5 X. Liu, H. Hu, P.D. Drummond, Phys. Rev. A78, 023601 
(2008); P. Kakashvili and C.J. Bolech, Phys. Rev. A79, 
041603(R) (2009); T.K. Koponen et al, New J. Phys. 10 
045014 (2008). 

6 M. J. Wolak, V. G. Rousseau, C. Miniatura, B. Gremaud, 



5 



R. T. Scalettar, G. G. Batrouni, larXiv: 1004.44991 Phys. 
Rev. A82, 013614 (2010). 

R. Blankenbecler, D.J. Scalapino, and R.L. Sugar, Phys. 
Rev. D24, 2278 (1981); S.R. White et ai, Phys. Rev. B40, 
506 (1989). 



V.G. Rousseau, Phys. Rev. E 77, 056705(2008); V.G. 
Rousseau, Phys. Rev. E 78, 056707(2008). 
G. G. Batrouni, M. J. Wolak, F. Hebert, V. G. Rousseau; 
Europhys. Lett. 86 47006 (2009). 



